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Abstract. This paper is concerned with diffuse- interface approximations of the Willmore flow. 
We first present numerical results of standard diffuse-interface models for colliding one dimen- 
sional interfaces. In such a scenario evolutions towards interfaces with corners can occur that do 
fS| not necessarily describe the adequate sharp-interface dynamics. 

We therefore propose and investigate alternative diffuse- interface approximations that lead 
to a different and more regular behavior if interfaces collide. These dynamics are derived from 
approximate energies that converge to the L^-lower-semicontinuous envelope of the Willmore 
O ^ energy, which is in general not true for the more standard Willmore approximation. 
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1. INTRODUCTION 



Diffuse-interface approximation of geometric evolution equations has a long history and is 
I— I widely used in numerical simulations. One advantage of the diffuse-interface approach is that 

^ usually an automatic treatment of topological changes in a reasonable manner is guaranteed. It 

is however not always clear which (generalized) sharp-interface evolution is in fact approximated. 
4^ We discuss here this issue in the case of the diffuse Willmore flow, in particular in situations 

^ where collisions between different interfaces occur and where these interfaces interact with each 

^ other. In applications such as image processing and computer vision, it is of great interest to 

I— I compute Willmore flow through such topological events. 

^ In the following we fix a nonempty open set C M^. Let A4 denote the class of open sets 

^ E C ft with r = dE n Q given by a finite union of embedded closed (n — l)-dimensional C^- 

^-H manifolds without boundary in ft. We associate to such F the inner unit normal field z/ : F ^ R^, 

the second fundamental form A with respect to z/, and the principal curvatures ^i, . . . , Hin-i with 
\^ respect to u. Finally we define the scalar mean curvature i7 = /^i + . . . + hin-i and the mean 

curvature vector H = Hu. 
Q The Willmore energy [38j is then defined as 

(N 

Tl mr):=;^ / H\x)dn^-\x). (1.1) 

The corresponding L -gradient flow is called Willmore flow. For an evolving family of sets 
^ (^(^))tG(o,T) ill with boundaries F(t) = dE{t) nri the velocity in direction of the inner normal 

field v{t) is given by 

v{t) = Ar(,)if(t) - \H{tf + H{t)\A{t)\\ (1.2) 

on F(t), where |^(t)p = + • • • + denotes the squared Frobenius norm of the second 

fundamental form A{t) and ^T{t) denotes the Laplace-Beltrami operator on F(t). 

In the case of two space dimensions the Willmore functional for curves and the Willmore flow 



are better known as Eulers elastica functional and evolution of elastic curves. In this case (1.2) 
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reduces to 

v{t) = Anit) + ^K{tf, (1.3) 

where Hi(t) denotes the curvature of T(t). The Willmore flow of a single curve in the plane exists 
for all times [23] and converges for fixed curve length to an elastica. 

A well-known and widely used diffuse-interface approximation of the Willmore energy is moti- 
vated by a conjecture of De Giorgi [18j and is in a modified form given by 

We{u) [ ^( -£Au+-W\u)y dC'', (1.4) 

where W' is the derivative of a suitable double- well potential W and is a smooth function on R^. 
The corresponding formal approximation of the Willmore flow is then given by the L^(fi)-gradient 
flow of Ws, 

edtu = -A(-5A^x+ ^H^'H) + ^H^''H(-5Ai/+ Jh^'H), (1.5) 

complemented by suitable boundary conditions for u on dVt and an initial condition for i/, in fi. 

If contact and collision of (sharp) interfaces are possible it is a priori not clear how to continue 
the Willmore flow. This situation already occurs for the evolution of several curves in the plane. In 
many applications interactions of different curves should be considered and treating the evolution 
of each curve separately might not be appropriate. To account for touching and colliding curves 
in such situations a suitable generalization of the evolution beyond the smooth embedded case 
is required. One possible extension is a gradient dynamic with respect to a suitably relaxed 
Willmore functional for general sets. A natural candidate for such relaxation is the L^(fi)-lower- 
semicontinuous envelope 

y^{E) := inf (liminf>V(a£;fc) : E^^^ E in L^{n),Ek E M for ah k G n|. (1.6) 

It is however difficult to characterize and numerically approximate the corresponding gradient 
flow. Instead we consider here diffuse-interface approximations that naturally allow for collision 
of (diffuse) interfaces and exist globally in time. 



Our first observation is that the usual diffuse approximation (1.4) does not approximate any 



gradient flow with respect to the generalized elastic energy (1.6). In fact, careful numerical 
simulations show that the diffuse evolution in the plane can lead to transversal intersections of 
interfaces, which is known to have infinite energy with respect to the L^(R^)-lower-semicontinuous 
envelope of the elastica functional [8] . The occurrence of intersections is related to the existence 
of saddle solutions of the Allen-Cahn equation, as we will explain below. This example also 
demonstrates that already on the level of energies a discrepancy between the diffuse elastica 



energy and the relaxation (1.6) of the sharp-interface energy occurs. In particular. We does not 
Gamma-converge to VV [32j. 

We will next discuss two diffuse-interface approximations of the Willmore flow that Gamma- 
converge to VV. The first one was proposed by Bellettini [5j. As an alternative we introduce 
a modification of the usual diffuse Willmore functional that adds a penalty term to W^. This 
penalty term is very small as long as no collisions of interfaces occur, allows for the touching of 
diffuse interfaces, but prevents transversal intersections. For both approaches we consider the 
L^-gradient flows and present numerical simulations. 
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2. Diffuse approximation of the Willmore functional and diffuse Willmore 

FLOWS 



Diffuse-interface approximations of the Willmore energy ( 1.1 ) are usually based on a Ginzburg- 
Landau free energy for a phase-field variable u 

ne{u) := ^ + ^~^W{u)^ dx, (2.1) 

where W denotes a suitable double- well potential that we choose in the following as W{u) = 
18v?{l — u^. For a smooth phase field : ^ R let us further define the L^-gradient of 7/^, 

ou 

the diffuse normal vector 



else, 



and the level set mean curvature 



v{x) := V-z^(x). (2.2) 

For phase fields u with 'moderate' energy T-Le{^) the function u will look like a smoothened 
indicator function that is close to the values 0, 1 in a large part of the domain and possibly 
with thin transition layers. The width of these diffuse interfaces is proportional to £ > 0. 
The energy 1-Le is a diffuse-interface counterpart of the surface area functional. This statement 
was made precise by Modica and Mortola, who proved the Gamma-convergence with respect 
to of T-L^ to the perimeter functional [3ll The L^-gradient w oi describes a kind 

of diffuse mean curvature and motivates the definition of the diffuse Willmore functional ( |1.4| ). 
The approximation We has been studied intensively and is widely used in numerical simulations 
[27t O [221 EH]- For space dimension n = 2,3 it is known [33j that (for uniformly bounded He) 
the functionals We Gamma-converge towards the Willmore functional in limit points E G ft with 
(7^-boundary in fi. The Gamma-convergence is not true in general limit points E C M^. In fact 
[T6j showed the existence of a smooth function : ^ (— 1, 1) (where ±1 are the zeros of the 
double well potential) with the following properties: 7i is a saddle solution of the Allen-Cahn 
equation 

-Au + W\u) = inR^ 

u = holds on the coordinate axes {(xi,X2) : xiX2 = 0}, and u is positive inside the first 
and third quadrant of R^ and negative inside the second and fourth quadrant. Moreover, there 
exists A: > such that for \y\ > k the saddle solution u is exponentially close to ±1 and 
\Vu\ is exponentially small. In particular, the rescaled functions (x) := u{^x) 

have on compact subsets of R^ uniformly bounded energy 7^^, converge in 

L\ocO^^) to the ±1- 

characteristic function of the set E := {(xi, X2) : ^1X2 > 0}, and finally satisfy We{ue) = for all 
£ > 0. On the other hand, we have by [8J that W{E) = 00 for the LjQ^(R^)-lower-semicontinuous 
envelope of the Willmore functional. Therefore 

oc = W(E) > liminf H^^(^x^) = 

which shows that We cannot Gamma-converge to W. 

For a tighter convergence of approximations some control on the Willmore energy of level sets 
of the phase fields is necessary. Bellettini [5j proposed such type of approximations for general 
geometric functionals. The corresponding approximation of the Willmore functional is given by 
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the square integral of the mean curvature of the level sets of the phase field u integrated with 
respect to the diffuse area density, 

mu):=l [ (v^)V^|V.xp + s-M^))dx. (2.3) 



In the particular case of (2.3) Bellettinis results imply that (for uniformly bounded diffuse surface 
area l-L^) the functionals in fact Gamma-converge with respect to L^{VC) to VV. 

For n = 2 an alternative approximation of the elastica functional has been investigated by 
Mugnai [32j. He uses an approximation of the square integral of the second fundamental form, 

dx. (2.4) 



\/u\ \\/u\ 

and obtains for n — 2 the L^(fi)-Gamma-convergence of VV^ to VV, again under a uniform bound 
on the diffuse surface area. 

3. A NEW DIFFUSE-INTERFACE APPROXIMATION OF THE WiLLMORE FUNCTIONAL 

We propose here a modification of the 'standard' approximation of the Willmore functional 
by an additional penalty term. This has some advantages in numerical simulations as we discuss 
below. The example of saddle solutions for the Allen-Cahn equation and a comparison with 
We-^ySJe reveals that the standard approximation works well as long as phase fields u behave like 
the optimal profile q for the one-dimensional transition from to 1 given by 

-q ^W{(i) = 0, g(0) = 0, lim q{r) = 0, lim g(r) = 1. (3.1) 



-oo 



Formal asymptotic expandions often use that u^{x) ^ ^(f ), where d denotes the signed distance 
functions to a limit hypersurface. This property can be formalized as vanishing of the discrepancy 

Q:^^-\Vu\^-s-'Wiu) 

that measures deviation from equi-distribution in the diffuse surface energy. In the limit £ ^ 
this quantity vanishes in L^(Q) for sequences of phase fields that have uniformly bounded diffuse 
surface energy and diffuse Willmore energy We [33j . This weak control does however not exclude 
formation of transversal intersections in the limit. Rewriting the diffuse mean curvature in terms 
of the level set mean curvature := V • as 

w -eAu + -W^u) -elVulv - VC^ • t^^^^ on {Vu ^ 0} 

we see that the diffuse mean curvature controls the level set mean curvature if and only if the 
normal projection of the gradient of the discrepancy is small. 
This motivates to introduce a penalty term of the form 

Mu) := ^ l^(w+ [e\Vu\ V2Wiu)) ' dx, (3.2) 

where < a < 1. If g(f ) we find that Aeiu) = 0{e) remains small. The energy however 
becomes large if u deviates from the optimal profile structure. We then define the modified diffuse 
Willmore energy 

:F,{u) := We{u)^Ae{u). (3.3) 
Using Bellettinis result [5] we can prove the Gamma-convergence of the modified Willmore energy 
J^e to VV. In the following we extend the functionals T-L^ and J^^ to L^(fi) by setting them to +oo 
on Li(J^) \C2(J^). 

Theorem 3.1. Let a > 0. Then the functional Gamma-converges with respect to L^{ft) to 
W in the following sense: 
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(1) Let {u^)^yo be a sequence of smooth phase fields : Q ^ with 



sup 7^5(^X5) < oc (3.4) 

e>0 



and ^ u in L^{Q). Then u G BV{Q.] {0, 1}); and 



W{E) < liminf J'^(^x^), (3.5) 



where E = {u = 1}. 



(2) Let E d he given with W{E) < oc. Then there exists a sequence {u^)^^^ of smooth 
phase fields in such that Xe in L^{Vt) and 

W{E) > limsnpTeiue). (3.6) 



£^0 



In the situation of ([!]) for a = we still obtain 

1 



W{E) < limmiT,{u,) (3.7) 

2 e^O 



and in particular }V{E) < 00 if the right-hand side is finite. 

Proof. (1) By Young's inequality we deduce for any K > the following estimate. 

/ - \ ^ 

w^ + K (w + (^e\Vu,\^/2W^y v) 

= (1 + K)w^ + 2Kw [e\Vu,\^/2W{ue)^ ' v + Ke\Vue\y^2W^e 



>e-^\Vu,\^2W{u,y. (3. 



1 + K 
This implies 



C ^2W{s) [ v^{x)dn''-\x)ds, (3.9) 



e 

> 



2 + 26- 

where we have used the co-area formula in the last line. 

For {us)e>o with ^ u in L^{Q) we first have by the Modica-Mortola Theorem that 
u E BV{n; {0, 1}) with 

/ \Vu\ < liminfHJi/^). 
Jn 

After passing to a subsequence (ek)keN we may assume 

lim J^ek{ue^) = liminf J'^('^e). 

k^oo £—^0 



By (3.9) we now can follow the proof of [5^ Thm. 4.2]. First one obtains a subsequence 
k ^ 00 (not relabeled) and a set / C (0, 1) with full measure such that for any 5 G / 

{usk = 5} = d{ue^ > 5}, 
{u,, =s}n {Vu,, = 0} = 0, 
^{uej^ys} ^{u>s} = Af^; as A: ^ 00, 
where E = {u = 1}. Moreover by the definition of W we have 



W{E) < liminf>V(A'|^^ >a) = liminf / dW'^ 
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for any s E I. By (3.9) and Fatou's Lemma for a > we eventually obtain 



liminf >W{u) [ ^y2W{s)ds = W{u). 

Jo 



In the case a — the same argument shows (3.7). 
(2) Let first E C ^ have smooth boundary. We follow the standard construction of a recovery 
sequence and consider the one-dimensional optimal profile q from ( |3.1| ) and the signed 
distance function d from dE. We then set := q{^) in {\d\ < 6} where 5 > is suitably 
small such that the projection on dE is smooth on {\d\ < 26}. In {\d\ > 26} we set to 
1 in and in ^}\ E. In {6 < \d\ < 26} we choose to smoothly interpolate in such 
a way that yV{ue) and \Vus\ are exponentially small in {\d\ > 5}, see for example [201 
Section 4]. Then the Willmore energy W£{ue) is known to converge to W{E). For the 
additional part in the energy we compute in the set {\d\ < 6} 



w + (e\Vu\^/2W^Y V = -q\^)Ad+{^J2W{q{^))q\^))'2v = 0, 



since in < 5} we have y^2W{q) = q^ and v = Ad as level sets of correspond in that 
region to level sets of d. 

In the region {\d\ > 26} we have e\\/u\^/2W{u) = 0. Finally, for a carefully chosen 
interpolate in the construction of u^, in {6 < \d\ < 26} we have that w and e\\/u\y^2W{u) 
are exponentially small. Furthermore, Ad = v is controlled in terms of the principal 
curvatures of dE. This shows the approximation property for C^-boundaries dE. 

To deal with the general case we only need to show that VV = W on sets with C^- 
boundaries. This condition is in fact satisfied by [35] and [29] . 

□ 



4. Numerical simulations for the standard diffuse Willmore flow 

Numerical investigations of the Willmore flow and related dynamics mainly build on parametric 
(sharp-interface) approaches and implicit treatments by level set and phase-field methods. Para- 
metric approaches for the Willmore flow have been proposed in [21 H] for curves and in [HU [191 E] 
for curves and surfaces. Generalized Helfrich-type flows for single- and multicomponent vesicles 
have been studied with sharp-interface methods in [13j and [25j, respectively. Level set methods 
have been applied first in [21] to the Willmore flow. For a comparison of level set and sharp- 
interface approaches we refer to [10]. Phase field approximations for the Willmore flow have 
been numerically investigated in [22]. For diffuse-interface approximations of Helfrich-type flows, 
we refer to [T2l [23l [Mj . Coupled Helfrich- and Cahn-Hilliard-type flows have been numerically 
treated with phase-field models in [371 128]. 



In this section, we consider the standard diffuse Willmore flow (1.4) and focus on situations 



where (diffuse) interfaces collide. We present numerical results for both a finite element dis- 



cretization and a finite difference scheme of the diffuse-interface flow (1.5). 



4.1. Finite element approximation. The discretization is implemented in the FEM toolbox 
AMDiS [36j. First, we use linear elements in space and a semi-implicit time discretization with 



a linearization of nonlinearities. Furthermore, we use a uniform grid and discretize (1.5) as a 
coupled system of two second order PDE's for the discrete solutions and and solve the 
resulting linear system with a direct solver (UMFPACK, [17j). Thereby, we consider a domain 
= ( — 1, 1)^ C and assume periodic solutions at the boundary dQ. Moreover, we use a simple 
adaptive strategy in time, where time steps Atm E [10~^,5 • 10~^] are inversely proportional to 
the maximum of the discrete time derivative of the phase-field variable. For the results presented 
in this section we have used £ = 0.1. 
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4.1.1. Symmetric initial condition. As initial condition, we have chosen a phase- field function 
Uhi'^O) : O ^ R having nine symmetrically distributed circular levelsets {uh{-, 0) = 1/2} (Fig. [l| 
left) with equal radii 0.1. In Fig. [l| one can see the contour plots of Uh at different times. Thereby 
discs start to grow until the interfaces begin to "feel" each other. Then the interface forms sharp 
corners. Our interpretation of such behavior is that the diffuse approximations converge to the 
saddle solution of the Allen-Cahn equation, discussed in Section [2j As the diffuse mean curvature 
for saddle solutions vanishes, such corners carry very little diffuse Willmore energy. 
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Figure 1. Evolution of "standard" diffuse- interface Willmore flow ( |1.5| ): Discrete phase-field 
Uh for different times t = 0, t ^ 0.0019, t ^ 0.0024, t ^ 0.0037. 



4.1.2. Non-symmetric initial condition. In Fig. [2j similar phenomena can be observed for non 
symmetric initial conditions. The discrete Willmore energy is plotted in Fig. [3} For large times 
t this energy 

is close to 0. 
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Figure 2. Evolution of "standard" diffuse- interface Willmore flow (1.5): Discrete phase-fleld 
Uh for different times t = 0, t 0.0039, t ^ 0.1038, t ^ 0.3338. 



4.1.3. Two circles initial condition. As a further example, we consider an initial condition with 
two circles with radii 0.2 and 0.3. Contour plots of at different times are shown in Fig. [4| 
where after collision of interfaces a transversal intersection appears. 

4.2. Finite differences implementation. In this subsection, we implement the standard dif- 
fuse Willmore flow using finite differences, and observe the same unexpected behavior in 2D as 
in the previous subsections: the colliding interfaces form a "cross", even though this should be 
precluded by the sharp-interface energy. Therefore, it seems unlikely that this surprising behavior 
is due to a numerical artifact: it appears to be an intrinsic feature of the standard approximation 



(1.4). However, we also provide simulations in three dimensions that suggest that this phenome- 
non might be two dimensional only: The standard diffuse-interface model inspired by De Giorgi's 
conjecture appears to lead to a topological change (mergers) when two typical surfaces collide 
in 3D. This is in keeping with some of the numerical experiments carried out by Du et. al. in 
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diffuse Willmore energy 




Figure 3. Evolution of "standard" diffuse- interface Willmore flow (1.5): Diffuse Willmore 
energy e(t) = yVe{uh{-,t)) versus time t. 



o 



o 




Figure 4. Evolution of "standard" diffuse-interface Willmore flow (1.5): Discrete phase- field 
Uh for different times t = 0, t ^ 0.003, t ^ 0.0045, t ^ 0.0200. 



j. Whether Gamma-convergent numerical approximations, such as the one due to Behettini 
[6] or the new one presented in this paper in Section [sj also lead to a topological change in 3D 
under these circumstances will be investigated subsequently in Section [5j See also Section [7] for 
a discussion of why topological changes are in fact more likely in 3D than in 2D for gradient 
descent of the relaxation of Willmore energy. 

In 2D, the results presented in Figure [5] are for the diffuse-interface energy 

where the second term is included to ensure a uniform bound of the diffuse surface area (note 
however that for evolutions with well-behaved initial conditions this is often automatically satis- 
fied, as for example in the simulations above). 
Our scheme for its gradient flow is: 

= -Ahw''+\^W'\u'') + -f\w'' (4.1) 



6t 



where 



w 



(4.2) 



and Ah is the standard centered differences discretization of the Laplacian on a uniform grid. 
This is an explicit time stepping scheme, the stability (CFL) condition (upper bound on the time 
step size 6t) for which scales as (5x)^ as 6x^0^. 

Alternatively, we have the following semi-implicit version: 

^,72+1 



St 



w 



(4.3) 
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where is again as in (4.2). At each time step 



,n+l 



is solved for via the discrete Fourier 



transform. This scheme appears to be stable for much larger time step sizes than (4.1). In the 
interest of minimizing the possibility of numerical artifacts, we refrain from further attempts to 
improve the computational efficiency of the schemes used here, even though there is no shortage 
of classical techniques for doing so. 

Numerical simulations using scheme (4.3) are shown in Figure [sj The computational domain 
was [0, and the diffuse-interface parameter was chosen to be s = 0.005. The spatial res- 

olution was 200 X 200. The parameter 7 was taken to be |. The results testify to the same 
surprising qualitative behavior as in simulations of Sections 4.1.1 and 4.1.3 Thus, this appear- 
ance of crossings when interfaces collide appears to be a robust, inherent feature of the standard 
diffuse-interface approximation. 









00 


00 



Figure 5. Gradient flow for the standard diffuse- interface approximation of WiUmore energy, 
using the finite differences scheme [4.31 



We now turn to some 3D simulations, again with the standard diffuse-interface approximation 



(1.5). The natural analogue of the two disks initial condition in three dimensions is two disjoint 
spheres. However, unlike disks in 2D, spheres in 3D are stationary under WiUmore flow, and would 
in fact shrink to naught in the presence of even the slightest additional penalty on perimeter (i.e. 
when 7 > 0). We therefore add an expansionary bulk energy term: 



yVe{u) + ^T-Ls{u) — a / udx. 



(4.4) 

case, adapt trivially to the 



with a > 0. Schemes (4.1) and (4.3), which correspond to the a 
a case. 

Figure [6] shows that the inclusion of the expansion term in 2D does not alter the formation of 
a cross (and failure to merge and become regular) for the two disk initial data. This simulation 
was carried out with the same parameters as before and a = 10. 






00 


00 


00 



Figure 6. WiUmore flow with perimeter penalty and volumetric expansion term, simulated 
using a finite differences discretization of the standard diffuse- interface approximation. 

Figure [7] shows simulations in 3D with two disjoint spheres of equal size as initial data. The 
parameters were e — 0.01, 7 = |, and a — 10. The spatial resolution was 100 x 100 x 100. 
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The volumetric term leads to expansion of the spheres, which eventually touch. Unlike the 2D 
situation, the two surfaces merge and become instantaneously regular. 



I 



Figure 7. 3D simulations with the standard diffuse- interface approximation of Willmore flow, 
together with a volumetric expansion term. 



5. Numerical simulations witr Bellettini's approximation 
In this section, we provide numerical experiments in the plane with Bellettini's Gamma- 



convergent approximation (2.3), and investigate what happens at topological changes. In the 
interest of isolating potential numerical artifacts from inherent behavior of Willmore flow that 
results from this approximation, we keep the numerical scheme as simple as possible: It is a 
straight-forward finite differences discretization, with explicit time stepping. Unsurprisingly, com- 
putations with this scheme are very slow, owing to the extremely stringent stability restriction on 
the time step size. However, they appear to be also very robust. In particular, a discrete form of 
the energy is observed to decrease at every time step. Due to the delicate nature of the question 
(topological changes in a fourth order gradient flow for a curvature dependent functional!), we 
give full details of the implementation. 



We will work with the following regularized version of (2.3): 



-Mi 



2 



o / { === \Vu\' + -W{u) dx. (5.1 



^\Vu\^ + 5 V2 



e 



Here 7 > is added to ensure a uniform bound for T-Le{u^), necessary for Gamma-convergence 
to the L^-relaxation of the Willmore functional. For simplicity, our exposition is restricted to 
below; extension to arbitrary dimensions is straight forward. To begin with, the gradient flow 



for ( |5.1[ ) leads to the following evolution: 

{u^ + 5)dx{i^6he) - UxUydy{nshe 




where 



(|Vix|2 + 5)2 ) 
{ul + 5)dy{hishe) - UxUydx{l^5he) \ (5.2) 

(|Vix|2 + 5)i I 
{nl + ^)Vu)-^y, + ^)W\u). 



^^-V-f^^^l (5.3) 



and denotes the diffuse surface area density, he{u) = ||Vi/p + s~^W{u). 

Working on a uniform spatial grid with periodic boundary conditions, let and D~ denote 
the standard forward and backward difference quotients in the direction of their subscript. Let 



APPROXIMATIONS OF WILLMORE FLOW 



11 



u'^ denote the solution at the n-th time step. Let's start with the discretization of the diffuse 
surface area density h^(u): 

M{u) = I [{Dtuf + {D-uf + {D+uf + {D-uf) + hv{u). (5.4) 

The denominator in the curvature term k.^ can be discretized in any one of the fohowing four 
ways (we'h use ah): 



Du^^^ = ^{dUY + [Dfuf + 5 (5.5) 

where the first and second zh in the superscript of Du refer to the signs of difference quotients in 
the X and the y directions, respectively. The superscript will be dropped for convenience below, 
whenever it is just ±, ±. The choice of sign for the difference quotient in each coordinate direction 
in (5.5) determines the signs of all subsequent difference quotients, as indicated with zb or =F signs 
below. Approximation to the curvature term can be obtained as: 

(5.0) 

where the superscripts of K indicate the signs for the difference quotients in the x and y directions, 
respectively. 



Next, define 



[Dt{K{u)M{u))] [{D^uf + S]\ ( [DtiKM)] {Dtu){D^u) 



A^,Hu) := I '"^ I - 




± 



u) 



{Duf 

[D^{K{u)M{u))] {Dtu){D. 

{Duf ) ■ 

(5.7) 



where superscripts of Ai and A2 indicate once again the chosen signs for the difference quotients 
in the x and y directions, in that order. Let 



(5.8) 



A\{u) := Af'^iu) + At'-{u) + ^^'+(n) + ^^'"(n), 
A^u) := At'+iu) + At'-iu) + A-'+iu) + A-'-{u). 
Define also the discrete squared curvature: 

kUu) := ] {{K+'+{u)f + {K+'-iu)f + {K-'+{u)f + {R-'-iu))') . (5.9) 



Let 



Mu) ■.^-l{D-{iK^,+j)D^u^)+D;{{K^,+j)D^u^) + 

Dt{{K^,+^)D-u^)+D+{{K^,+^)D;u^)} (5.10) 



Finally, our update scheme is: 



^Ai{un + A2{un + -A3{un). (5.11) 

With time step size 5t > chosen small enough compared to the spatial grid size, this scheme is 
guaranteed to decrease the following discrete form of the energy 

Y,{K'.Kj)+7)M{uIj) (5.12) 
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as can be easily verified by differentiating (5.12) with respect to time, and summing by parts a few 
times. Of course, it must be mentioned that the convergence of energy (5.12) to, say, (5.1) as the 
grid size and the regularization parameter 6 are appropriately sent to has not been established 
rigorously; we merely give some numerical evidence. 

In 2D, the computational domain was [0, -^]'^ with a spatial resolution of 100 x 100. Periodic 

boundary conditions were used. The parameters were chosen to be 5 = ~ \^ ^ ~ 0.01. 

As an initial guess, vP was taken to be the characteristic function of the union of two disjoint 
disks. During the evolution, the disks are observed to initially expand as disks, as expected. 
However, once they are within a small distance (related to the diffuse-interface thickness e) of 
each other, they are unable to get any closer. Topological changes appear to be precluded; in 
particular, neither a merger subsequently leading to a smooth evolution occurs (as was previously 
seen in various numerical implementations of Willmore flow via implicit representations), nor a 
corner forms as with the De Giorgi approximation. Instead, the curves continue to expand, but 
are no longer circles. 




Figure 8. Evolution of disks under gradient flow (in the bulk) for Bellettini's approximation 
to the Willmore energy. Topological changes appear to be precluded: the disks do not merge. 
After getting within a small distance of each other (related to diffuse-interface thickness), they 
continue their expansion, but are no longer circles. 



Next, we explore what happens in 3D. Once again, as in Section [42| we include a volumetric 
expansion term so that spheres would grow. Figure |9] shows computations with two parallel 
cylinders as initial condition, which in fact corresponds to a 2D computation; the difference from 
the experiment of Figure |8] is the inclusion of the volumetric expansion term. We see that even 
in the presence of the additional driving force bringing the interfaces into collision, topological 
change is precluded, just as in the 2D experiment of Figure [8| 



On the other hand. Figure [TO] shows simulations with two disjoint spheres as the initial con- 
dition. The computational domain is [0, -^=\^ and the spatial resolution is 50 x 50 x 50. The 

parameters were chosen to be 5 = 0.02, 7 = |, and 5 — 0.01. In this case, the spheres merge 
once they come into contact (i.e. "feel" each other due to the diffuse-interface thickness). See 
the remark in Section [7| for an explanation of why this topological change is not precluded. 

6. Modified diffuse-interface Willmore flow 



In the following, we consider the modified diffuse approximation of the Willmore energy (3.3) 
and the corresponding L^-gradient flow. 

This modification by an additional 'penalty' term offers some extra flexibility compared to 



alternative energies such as (2.3). As long as solutions are close to the optimal-profile construction 
the correction is negligible. One could therefore 'switch off' the additional term as long as 
Ae is small, or only use the extra forcing in regions where Ae is large. This helps to reduce 
computational costs, as the remaining standard term is much easier to deal with. If solutions 
start to deviate from the optimal profile, we only use the property that the extra term blows up. 
We therefore need much less accuracy in computing this term and can choose different numerical 
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Figure 9. Evolution of cylinders under gradient flow for Bellettini's approximation of Will- 
more energy, with a bulk energy term encouraging expansion added. The evolution is essentially 
2D; the only real difference from the simulation of Figure [Slis the inclusion of the bulk energy 
term. Just like in the no bulk energy term case of Figure[8] the cylinders cannot merge even 
though they come into contact, despite the additional driving force that slams interfaces into 
each other. 










Figure 10. Evolution of sph eres under gradient flow for Bellettini's approximation of 
Willmore energy, with a bulk energy term encouraging expansion added. The spheres expand 
and touch, and unlike in the 2D experiment with disks shown in Figure [S] or the 3D experiment 
with cylinders shown in Figure |9] the topological change takes place: the spheres merge, and the 
surface becomes instantaneously regular. 



relaxation or discretization parameters than for the term We in the total energy J^^. Here however 
we only aim at a proof of concept and do not exploit such possibilities. 

6.1. Modified Willmore energy. In order obtain a simpler variational derivative than we 
would obtain from (3.2), we assume \e\/u\ ^ ^/2W(u) and introduce an additional energy 

Asju) ( -sAu + e-'W\ul + ^2H^V • ) dx, (6.1) 

= :w ^ V ' 

= :v 

where a < 1. Whereas the analysis in Section [3] does not cover this case, our numerical simulations 
show that this choice works equally well. In the following we therefore will use the diffuse Willmore 
energy 

and will analyze numerically this functional and t he c orresponding L^-gradient flow. We expect 



that a Gamma-convergence result as in Theorem 3.1 also holds for J^^, but the analysis in this 



case is much more difficult and will be subject to future investigations. 



6.2. Computation of the additional energy. We compute the additional energy 

a{t) :=Ae{uh{',t)) 
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with a = during a simulation similar to the one in Fig. [Tj The only difference is that we have 
used 4th order finite elements. Furthermore, we have regularized 

1 1 



with 5 = 0.1 in order not to divide by zero. Again we use a uniform grid with 65^ = 4225 vertices 



leading to 66049 degrees of freedom for each unknown. In Fig. 11 (left), one can see the energy 
decrease of e{t) = W£{uh{-^t)) during time with values very close to zero for large times. Fig. 
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(right) displays the expected blow up of the additional energy a{t) = A£{uh{-,t)) when the 
interfaces begin to "feel" each other. This behavior serves as a motivation to study the modified 
flow of the energy We + Ae in the following. 



diffuse Willmore energy 



additional energy 




2500 



2000 



1500 



500 





a(t) 


J 





0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 
t 



0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 
t 



Figure 11. Evolution of "standard" diffuse- interface Willmore flow (1.5): Willmore energy 



e(t) = yVe{uh{-,t)) (left), additional energy a{t) = Ae{uh{-,t)) (right) versus time t. 



6.3. Modified fiow. Here, we consider the L^-gradient flow 

du 



^9tu = ^- (6.2) 



of the energy W£{u) + A£{u). We introduce 

V := ^2W{u)V . ^ = V • (^V2W{u)^^ - {^/2W{^y\Vu\ 
and obtain the variational derivative 

^ = ^ f - sA{v + ^) + s-^W'{u){v + w)+ + ^) (6,3) 

+ V • {B{u, Vu){\/v + Vw))^ 

where 

\Vu\ \ |Vi/| |Vix| 

Thereby, / = denotes the unit matrix. 

In order to numerically treat the modified flow (6.2), we split the time interval [0, T] by discrete 
time instants = to < ti < • • • < = 7", from which one gets the time steps Atj^ := t^+i — t^, 
m = 0, 1,...,M — 1. Moreover, we apply the following operator splitting type Ansatz for time- 
discrete functions u^^\ v^^\ w^'^^ at time instant tm'. 
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(1) compute v^*") via 



I |Vn(™)| 
(2) solve for u^"'+^^ and in 

^(m+l) _ ^(m) 



Ml 



At, 



(6.4) 



(6.5) 



1 



rl+a 



£A(t;('") + ?/;(™+i)) + e-iW'(n('"))(t;(™) + w^"'+^'^] 



v/2W(nM) 
+ V • (S(n('"), Vn('"))(Vt;("^) + Vu'("^+i)))) 
y^(m+i) ^ _gAu('"+^) + £-^VF"(u('"))u("*+^) + e-^W'iu^"^^) - («("*) )«("*). (6.6) 

Thereby, for 5, 6w > 0, we regularize 



1 



1 



1 1 



v/2W(nM) ^2VF(nM) + 



for all terms of these forms appearing in ( |6.4[ )-( |6^ . In addition to the previous parameters, we 
have used 5w — 0.01. 

As a first test, we observe that the modified flow yields a reasonable approximation of Willmore 



flow in the case of a growing circle. In Fig. [T2j we compare numerical results with the analytic 
expression for a circle growing according to Willmore flow. Thereby, we use an initial radius 
i?(0) = 0.1 and plot the radius R{t) versus time. 



radius R(t) 




0.02 



Figure 12. Evolution of modiOed flow JO): Radius of growing circle versus time: Analytic 



expression and results of modified diffuse-interface flow. (6.2) 



In Fig. [T3| we see numerical results for this flow showing the phase-field variable u for a 
symmetric initial condition with nine circles with equal radii 0.1 as in Fig. [l} The final picture 
shows the nearly stationary solution. Fig. 14 shows the energy decrease for this example. 

In Fig. [15] the nonsymmetric initial condition from Fig. [2] has been used. Again, the nearly 
stationary discrete solution at time t ^ 0.016 shows that self intersections are prohibited for this 
modified fiow. 
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Figure 13. Evolution of modified flow (6.2): Discrete phase-field Uh for different times t = 0, 
t ^ 0.0012, t ^ 0.0024, t ^ 0.0045. 
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Figure 14. Evolution of modified flow (|6^: Total energy f{t) = We{uh{-,t)) + Ae{uh{-,t)) 
versus time t. 
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Figure 15. Evolution of modified flow (6.2): Discrete phase-field Uh for different times t 
t ^ 0.0010, t ^ 0.0030, t ^ 0.0160. 



0, 



In Fig. 16 the initial condition with two circles is not symmetric as in Fig. |4j The circles grow 
until the interfaces "feel" each other. The modified energy prevents the level sets from forming 
self intersections. The contour plots in Fig. 16 are taken at similar times as in Fig. |4j 

In the gradient flow simulations for Bellettinis energy and the modified energy J^^ we have often 
observed that circles upon collision keep touching and evolve to an ellipse type shape. We expect 
that such shapes represent suitable elastica. In a final example, we therefore compare (nearly) 
stationary states in our numerical simulations with graphs of minimal elastic energy [26j that 
present possible optimal configurations. In [26j, Linner and Jerome could prove the existence of 
a unique graph of minimal elastic energy among graphs that can be parameterized by y{t)) 
such that 

(x(0),y(0)) = (0,0) and {x\t),y\t)) = L(cos0(t), sin0(t)). 
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Figure 16. Evolution of modified flow (6.2): Discrete phase-field Uh for different times t = 0, 
t ^ 0.003, t ^ 0.0049, t ^ 0.0199. 




analytic minimizer [26] 
t = 
t - 2.2462 



>^ 




shifted analytic minimizer [26] 
t = 
t « 2.2462 



Figure 17. Evolution of modified flow ([6^: Level curve {uh = 1/2} at times t = and 
t ~ 2.2462 and analytic minimizer from [26] (left), level curve {uh = 1/2} at times t — and 
t ~ 2.2462 and analytic minimizer from [26] shifted by ^(1) in — ^/-direction (right). 



wfiere 9{0) = 0, 9{1) = | and (6>, L) G VF-^'^((0, 1)) x R+. Moreover, Linner and Jerome were able 
to provide an explicit formula for this minimizer. For the following comparison, we have chosen 
the rectangular domain Q = ( — 1.1, 1.1) x (—2.2, 2.2) with periodic boundary conditions such that 
we expect the level set {u^ = ^} to stop growing in the x-direction at x ±1. Furthermore, 
we chose an ellipse as initial condition. In Fig. 17 (left) we see level sets {uh = ^} at different 



times, where the nearly stationary discrete phaseneld Uh at time t ^ 2.2462 shows a reasonable 
approximation of the analytic minimizer from [26], which has been shifted by y{l) in — y-direction 



in Fig. 17 (right), indicating an even better approximation by the numerical result. 



7. Discussion 

We have analyzed the behavior of different diffuse approximations of the Willmore flow in 
situations where diffuse interfaces collide. Different scenarios emerge. In level set approximations 
the respective phases typically merge (not investigated here, but see the results presented in 
[21j) at the moment of collision. For the standard diffuse approximation on the other hand we 
have demonstrated that (at least in two dimensions) transversal intersections of diffuse interfaces 
are preferred. As an alternative we have introduced two new diffuse approximations with again 
different behavior: after collision the phases keep touching, away from the touching points the 
interfaces evolve to elastica. 
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Figure 18. Left: Curve 7fc and modification 7^. Right: Limit of the modified curves. 



Whereas all approximations converge to the same evolution as long as phases remain well 
separated, the kind of approximation determines the evolution past collision of interfaces. Here 
is a summary of these behaviors: 



Method 


Disks in 


Spheres in 


Standard 
Level set 
Bellettini 
New approx. 


"Cross" formation 
Merger 
No merger 
No merger 


Merger 
Merger 
Merger 

Not investigated (we expect mergers) 



Although there appears to be agreement in the case of spheres, it is reasonable to assume that the 
discrepancies in the two dimensional handling of topological changes would manifest themselves 
also in 3D with more general surfaces. 

The expected behavior through these topological changes, and therefore the choice of the 
approximation, may depend on the specific application. For example, in image processing appli- 
cations such as inpainting (see e.g. [llj), it is often desired that colliding interfaces merge in 2D. 
Indeed, in the important work [15j, the Willmore energy is utilized to drive interfaces towards 
collisions and mergers. In other applications, for instance in the segmentation of medical images 
(see e.g. [Ij), it may be desired to prevent merging of colliding interfaces. 

Our new approximation and that of Bellettini [5j agree in the handling of topological changes 
considered here for Willmore flow. Indeed, for these Gamma-convergent approximations, we argue 
that the observed limit flow describes the only reasonable evolution that is continuous with respect 
to the L^-topology and keeps decreasing the Willmore energy. In contrast, typical 'mergers' do 
not share this property. In fact, let us assume that we have a sequence of Jordan curves {'yk)keN 
that approximates the union of two touching balls 5 = {(x,y)GR^ : (x + l)+y^ < l}U{(x,y)G 
R^ : (x — 1) + 7/^ < 1} in L-*^- distance. Assume for simplicity that the 7/e are symmetric with 
respect to the x-axis, that the upper half is given as a nonnegative graph, and that the convergence 
IS m an -sense outside a region {(x,y) G R^ : |x| < 1 — 6}. We then can modify the curves 
by replacing the part below the x-axis by a circular arc that touches 7/^ in its intersection points 
with the X-axis (see Fig. 18). By construction the modified curves ^tre C^- Jordan curves that 
are piecewise C^. As the radius of the attached half-circles approaches 2 with k ^ 00 there exists 
a constant C > independent of k such that 



(7.1) 



On the other hand the sets enclosed by 7/e converge in L^-distance to (S'nR^) U (B(0, 2) n R?.) , 
where denote the upper and lower half-plane, respectively. This however is a set with exactly 
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one simple cusp point. By [71 Theorem 6.4] this imphes that 

hm yV{%) = oc. 



From (7.1) we deduce therefore that the elastica energy of blows up, too. 



In 3D the behavior of diffuse Willmore flows in situations where diffuse interfaces collide is 
different from the two dimensional case: two spheres that are initially disjoint but forced to come 
in contact via a volumetric expansion term are likely to merge under gradient descent for the 
relaxation of Willmore energy. Indeed, using e.g. catenoid "necks" , of the form 

y — a cosh 

V a 

with a and h are appropriately chosen, we can "connect" two spheres at the moment of contact 
with an arbitrarily small neck that is tangent to the spheres after slightly shifting them if necessary, 
while decreasing the energy. Note that the catenoid neck, regardless of its scale, has no elastica 
energy at all, since it happens to be a minimal surface and thus has vanishing mean curvature. 
To be more precise, the foregoing discussion shows that it is easy to construct a continuous in 
map from the interval [0, 1] into sets in that deforms the union of two spheres in contact to 
a surface topologically equivalent to the sphere, while decreasing the Willmore energy during 
the deformation. 
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